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r/. , Abstract. An efficient way of solving 2D stability problems in fluid 

O ■ mechanics is to use, after discretization of the equations that cast the 

* ^ ' problem in the form of a generalized eigenvalue problem, the incomplete 

^j, Arnoldi-Chebyshev method. This method preserves the banded siroc- 

co ■ ture sparsity of matrices of the algebraic eigenvalue problem and thus 

Oh' decreases memory use and CPU-time consumption. 

The errors that affect computed eigenvalues and eigenvectors are due to 

, ' the truncation in the discretization and to finite precision in the com- 

►> | putation of the discretized problem. In this paper we analyze those two 

^\ . errors and the interplay between them. We use as a test case the two- 

,-H ' dimensional eigenvalue problem yielded by the computation of inertial 

£N) , modes in a spherical shell. This problem contains many difficulties that 

make it a very good test case. It turns out that that single modes (espe- 
cially most-damped modes i.e. with high spatial frequency) can be very 
sensitive to round-off errors, even when apparently good spectral conver- 
gence is achieved. The influence of round-off errors is analyzed using the 
CZ3 ' spectral portrait technique and by comparison of double precision and 

extended precision computations. Through the analysis we give practical 
CZ) ' recipes to control the truncation and round-off errors on eigenvalues and 

^"\ eigenvectors. 
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1 Introduction 



The first step in studying the stability of the solutions of a nonlinear problem, 
is to solve the eigenvalue problem associated with infinitesimal pertubations 
which are superposed to the equilibrium state. Even if the equations of these 
perturbations are linear, solving the eigenvalue problem may be a formidable 
task. The difficulties arise in general when the variables do not separate. In 
such a case, the eigenvalue problem cannot be reduced to a set of smaller (one- 
dimensional) eigenvalue problems and one is left with a 2D or 3D problem. 



In most cases, after discretization of the equations, the temporal stability 
problem reduces to a generalized eigenvalue problem. A method to solve such 
problems is to use the QZ algorithm. Such an algorithm gives the full spec- 
trum of eigenvalues/eigenvectors, but the price to pay for obtaining this very 
rich information is very high in terms of memory requirement and CPU time 
consumption. Moreover the QZ algorithm does not preserve the sparsity of the 
matrices. On the other hand it is seldom needed to know the full spectrum: typi- 
cally one is interested in the few eigenvalues corresponding to the least stable or 
most unstable modes. For the foregoing reasons it is important to be able to solve 
the generalized eigenvalue problem with an iterative method which preserves the 
sparsity of the matrices and converges quickly and accurately to a small sub- 
set of the whole spectrum. Three types of iterative methods exist to solve such 
eigenproblems |8I9| . The first method is the subspace iterations method which is 
just a generalization of the well-known power method. A second method is the 
Jacobi-Davidson algorithm. The third one is to use Krylov based methods such 
as the Arnoldi method or the unsymmetric Lanczos method. Comparisons |11| 
between the subspace iterations method and the Krylov type ones tend to show 
that the second ones are more efficient when applied on large sparse matrices. 
We have chosen to use the Arnoldi method because it is easy to implement. Its 
backward stability is now well understood and it does not require any heuristics 
whereas numerical difficulties such as serious breakdowns can be encountered 
using the unsymmetric Lanczos method. 

In this paper we consider as a model problem the computation of the inertial 
modes of a rotating spherical shell. This problem contains many difficulties that 
make it a very good test case. A first difficulty is that the problem is essentially 
two-dimensional, because variables such as the radial distance and the polar an- 
gle cannot be separated. The size of the matrices thus grows very quickly with 
the resolution; for parameters of physical interest matrices are of order 10 5 or 
larger. The second difficulty is that the partial differential equations become 
of hyperbolic type and therefore yields, with boundary conditions, an ill-posed 
problem |14I13| . Solving this eigenvalue problem is therefore demanding numeri- 
cally. The third difficulty is that matrices are highly non-normal. The eigenvalue 
spectrum is thus very sensitive to machine precision and special tools must be 
used to analyze and control the round-off errors. Our analysis revealed some 
interesting aspects from the viewpoint of numerical precision. In particular we 
think that our results on the behaviour of round-off and spectral errors and their 
interplay are useful in many fields of physics where two-dimensional eigenvalue 
problems appear. 

We organized the paper as follows: we first describe the physics of our test 
problem and how we discretize it using spectral methods (section 2). We also 
briefly recall the principle of the incomplete Arnoldi-Chebyshev algorithm (sec- 
tion 3). We then discuss the role of spectral resolution (sect. 4) and presents 
our results about the behaviour of round-off errors (sect. 5); conclusions and 
outlooks follow. 



2 The test-problem 

2.1 Formulation 

We consider the problem of finding the modes of oscillation of a rigidly rotating 
fluid contained in a spherical shell; it is investigated for its astrophysical and 
geophysical applications (see |14I13| for a detailed discussion). Eigenmodes of 
this system are called inertial modes. 

The fluid is contained between two spheres of radii rjR and R (j) < 1) and 
rotates at an angular velocity fl around the z-axis. Choosing R as the length 
scale and (2J7) -1 as the time scale, the non-dimensional form of the equations 
for the equations governing perturbations are: 

J EAV x u — V x (e z x u) = AV x u , . 

where A is the eigenvalue (a non-dimensional frequency) and E = v/2QR 2 is the 
Ekman number. E is the non-dimensional measure of the kinematic viscosity v 
and is usually a small parameter (E < 10 -4 ). 

Equations Q are completed by boundary conditions on the velocity taken 
at r = r) and r = 1. We impose stress- free boundary conditions, namely that 

d (u s \ d 

u r = 



or \ r / or \ r / 



dr V r / dr 
on the boundaries (r, 0, ip are the usual spherical coordinates). 

2.2 Numerical method 

We discretize the preceding partial differential equations using spectral methods 
because of their efficiency at convergence |l(Jlci| . 

For obvious geometrical reasons, the angular part of the fields is expanded 
on spherical harmonics; hence, we set 

u = E E u U r ) R T + vL(r)Sf + w e m (r)TT, 

1=0 m=-£ 

where 

RT = Y e m (8, (f>)e r , ST = VYf 1 , Tf = V x Rf 

and where Y[ n (8,(j)) are normalized spherical harmonics (gradients in the defi- 
nitions of S" 1 and T™ are taken on sphere of unit radius) . 

Following some simple rules, given in J^l > the equation of vorticity (^i) may 
be projected rather easily on spherical harmonics. The radial functions u^r) 
and u^(r) then obey the following system 
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+A m {£+l)r~ 



dr 



r e+3 u i + 1 I = Xwi 



EA l + j^jMrul) - B m {iy- l ] 



dr 



■B m (l+\y- t -*±[r**«fi*\=\Mrn' w ) 



(2) 



where we have eliminated the v^s using V ■ u 
also been introduced : 



0. The following notations have 
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System (J2J is an infinite set of differential equations where the coupling between 
radial functions of indices £—l,£ and £+1 is due to the Coriolis force. Note that 
different m's are not coupled. 

For the discretization in radial coordinate we approximate the radial func- 
tions by truncated expansions of iV+1 Chebyshev polynomials. Thus, each of the 
functions may be represented either by its spectral components or by its values 
on the Gauss-Lobatto collocation nodes. We use the latter representation. In 
such case, differential operators d k /dr k are represented by full matrices of order 
(JV + 1). As system (J2J couples radial functions of indices t — 1, I and £ + 1, it 
yields a generalized eigenvalue problem with tridiagonal block matrices which 
we write symbolically: 
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For each value off, A I _ 1 e , Aq e and A\ £ , are (N+3) x (7V+3) matrices. They 
correspond to the discretization of the l.h.s. of the first equation in (J2J at the N+l 



Gauss Lobatto nodes, and of the boundary condition r 



dw'„ 
dr 



Wt 



imposed 



at the two radial boundaries. £ runs from m to L by steps of two when m is even 
and it runs from m + 1 to L by steps of two when m is odd. Similarly A l \ g , A l J- e 

and A 1 ^ are (N+5) x (N+5) matrices corresponding to the discretization of the 
l.h.s. of the second equation in @ at the Gauss Lobatto nodes, plus boundary 



conditions ui 



d A u* n 
dr 2 



2 rf"n, 
r dr 



at the two radial boundaries. 



3 The Incomplete Arnoldi-Chebyshev Method 

For efficiency reasons and memory requirements, the generalized eigenvalue prob- 
lem © should be solved using an iterative method because the matrices are large 
and sparse. As previously stated, a good method is the incomplete Arnoldi- 
Chebyshev method which we now briefly recall. 

Let K(u,A) = {u, Au, . . . , A m ~ l u] be the Krylov subspace built from the 
initial vector u, V m — {ui}i=i...m of size n x m be an orthonormal basis of this 
subspace. 

For applications to stability problems, one is mostly interested in the least- 
stable (or most unstable) eigenmodes which are associated with the generalized 
eigenvalues A with the greatest real part. Since these eigenvalues obviously do 
not belong to the outside part of the spectrum, we have to perform a spectral 
transformation. Let (n,y) be the solutions of 

([A]-a[B})- 1 [B]y = ^y. (4) 

Then, one easily shows that (A = <r+l//i, x = y). Thanks to this spectral 
transformation, the eigenvalues near the shift (the guess) a are sent to the out- 
side part of the spectrum and the Arnoldi method can now deliver the desired 
eigenpair very efficiently. The derived method can be summarized as follows 

Algorithm 1 Parameter: integers r (number of desired eigenpairs), m (number 
of Arnoldi steps), with r < m <C n, Arnoldi starting vector u and degree k of 
Chebyshev acceleration polynomial. 

1. Perform m steps of the Arnoldi method starting from u to compute V m 
and H m : 

{[A] - a[B])- l [B]V m = V m H m +v m+1 eJ n . 

2. Compute the eigenpairs (Hi,yi)i=i:m by applying the QR algorithm to H m : 

H m yi = [i l y l . 

3. If the stopping criterion is satisfied for the r wanted eigenvalues then go 
to step 6. 

4. Compute the parameters of the ellipse containing the m — r unwanted 
eigenvalues of H m and set zq = Y^i=i a iVm,yi where 

\\([A] - <j[B})- 1 [B})V m y l - HiV m yi\\ 



(\\[A]-*[B])-i[B]\\\\V myi \\ * 

5. Perform k steps of the Chebyshev acceleration starting from zq to obtain 
a better starting vector u for the Arnoldi method; go to step 1. 

6. Set {Xi = a + l//j^,xi= V m yi}i=v. r - 

This algorithm requires a matrix- vector product involving the matrix [B] and 
a linear solver to compute z-i solution of ([A] — a[B])z2 = z\. In our application, 
[A] and [B] are banded matrices, so that a band linear solver from LAPACK 
has been used. The internal dense eigensolver in step 2 has been taken from 
EISPACK. The interested reader is referred to |4l2lfil5| for more details. 



4 The role of spatial resolution 

We first study the convergence of eigenvalues and eigenvectors as a function of the 
resolution. The two relevant parameters are the degree of the largest Chebyshev 
polynomial (equal to the number of radial nodes N minus one) and the degree 
L of the last spherical harmonic. We deal only with axisymmetric m — modes 
and therefore drop the index m (no additional difficulty arises when m ^ 0). 

In the following we use the notation u> for the imaginary part of the eigenvalue 
(the frequency), and r for the real part. Thus A = r + icu. All the modes of our 
test-problem are stable, i.e. r < 0, and \t\ is the damping rate. 

We define the Chebyshev and Legendre spectra of the field u with spectral 
components u(£, n) in the following way: 



C(n) 



maxf \u(£, n)\ 
max^ n \u(£, n)\ 



£(£) 



max, 



\u(£,n)\ 



max^„ \u(£, n)\ 



Both spectra are filled because inertial modes display very fine structures 
(see J5 f° r typical spectra and eigenfunctions occurring at Ekman numbers as 
low as E = l(r 8 ). 

Here, we take a moderately small value of the Ekman number: E = 10~ 4 so 
that the full eigenvalue spectrum can be explored with an affordable resolution. 
As it may be expected, eigenvalues with smaller |r| require less resolution to 
converge than those with large damping rate (see figures |l(a)| and 1 (b) I . This 
is easily understood since eigenvectors with small |r| tend to have a smoother 
pattern, which is well approximated by a small number of spectral modes. 
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Fig. 1. Chebyshev (solid line) and Legendre (dashed line for ru and dotted line 
for w) spectra, (a): mode at E = 1CT 4 with w = 0.657976 and r = -0.00875. 
(b): mode at E = 10~ 4 with u = 0.654580 and r = -0.51. 



The convergence of eigenvalues as a function of spatial resolution goes to- 
gether with that of eigenvectors: unless all the scales present in the eigenvector 
are resolved, both the eigenvector and the eigenvalue are not well approximated. 
This observation allows us to give a simple rule to check the convergence of 
the eigenvalue. Let us define the ratio fi, between the lowest spherical harmon- 
ics coefficient and the largest one, and define gjy as the same ratio but for the 
Chebyshev expansion. 

min£(7) minC(n) 

■> L = 7ri\ 9n = -^-j—: (& 

maxl,({) maxC(n) 

These two ratios measure the truncation error in the spherical harmonic expan- 
sion (/i) and in the Chebyshev expansion (<?jv). We next define e as the absolute 
value of the difference between the computed eigenvalue and the converged one 
(i.e. obtained with a large resolution). 

In figure [21 we plot e as a function of //,. The number of Chebyshev polyno- 
mials was chosen large enough to resolve completely the radial dependence. We 
clearly see that e follows the law e oc f\ until a plateau is reached. The plateau 
appears at the largest resolutions and indicates that no better approximation 
to the eigenvalue can be obtained by increasing the resolution. It gives thus a 
measure of the round-off error of the computation. From the curves obtained 
for different eigenmodes, we see that the round-off error is a rapidly increasing 
function of the damping rate. 

In figure 13 we plot e as a function of the parameter gjq. Here the number of 
spherical harmonics was set large enough to fully resolve the angular dependence. 
The Chebyshev convergence appears to be governed by the law e (x g^. Here 
too, good convergence is obtained only for least-damped modes. We note also 
that the plateau values are very close to those of the preceding figure. 

5 The importance of round-off errors 

The foregoing results indicate that round-off errors play a major role in the 
accuracy of the numerical solution, especially for strongly damped modes. We 
shall now investigate this point more thoroughly. 

First, we stress the fact that good spectral convergence is not at all a guaran- 
tee against round-off errors. This point can be made very clear using the mode 



displayed on figure 1(b) for example. No doubts that for such a mode the spectral 
expansion has converged: there are 12 decades in the Chebyshev spectrum and 
16 decades in the Legendre spectrum; however the whole spectrum is subject 
to large round-off error at all wavenumbers. To illustrate this point we consider 
two different computations where we only change the value of the shift a (see 
equation Q) of the Arnoldi-Chebyshev algorithm: a — —0.51 + iO. 65458 in the 
first case and a = —0.51 + iO. 65558 in the second case. In both computations 
the Ekman number is E = 10~ 4 , L = 120 and N = 64; the Arnoldi-Chebyshev 
algorithm converges to the same eigenmode. The Chebyshev and Legendre spec- 



tra for the first case are those represented on figure 1 (b) the two spectra for the 
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Fig. 2. Error e of the computed eigen- 
value plotted as a function of the Leg- 
endre truncation error of the eigenvec- 
tor fi, (eq. (J5J ) . Different curves corre- 
spond to different eigenmodes. Thick 
line corresponds to the law e = j\. 
Note the horizontal plateau at large 
resolutions, due to round off errors. 



Fig. 3. Error e of the computed eigen- 
value plotted as a function of the 
Chebyshev truncation error of the 
eigenvector g^ (eq. (J5J). Thick line 
corresponds to the law e — gw- 



second case are similar. We plot in figure|Hthe relative difference of the spectral 
coefficients, defined as 



5C(n) = 



|C 2 (n)-Ci(n)| 

0.5(C 2 (n)+Ci(n)) : 



6C(n) 



|£ 2 (n)-A(n)| 
0.5(£a(n) + £i(n)) 



where subscript 1 (resp. 2) corresponds to first (resp. second) eigenvector. We see 
that the relative difference is spread almost uniformly throughout the wavenum- 
bers, until round-off error in the spectrum is reached, where necessarily the 
relative error grows to 0(1). This uniform spreading is not surprising; in £Q it is 
shown that for Chebyshev expansions the spectral round-off error of differential 
operators is distributed uniformly among wavenumbers. 

The round-off error may be investigated quite systematically by computing 
the spectral portrait of this eigenvalue problem. 

Spectral portraits and pseudospectra have recently attracted the attention 
as a tool of choice for investigating spectral properties of nonnormal matrices 
(see 171151161 1. It consists in the representation of the map 

z — > spp(z) = log 10 [||(,4 - zi?)- 1 || 2 (||A|| 2 + |*|||B|| 2 )] 

in a prescribed region of the complex plane. The contour lines of level e of the 
spectral portrait are the borders of the e-pseudospectrum of the matrix pair 
(A, B) : they enclose all the eigenvalues of the matrix pairs (A + AA, B + AB) 
and with ||AA|| 2 < e||A||a and IIABIU < e||5|| 2 . 
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Fig. 4. Relative difference of the spectral coefficients obtained with two com- 
putations where the only difference is a change of the Arnoldi-Chebyshev shift 



a of the Arnoldi-Chebyshev algorithm. The mode is that of figure 1(b) Dashed 
line: Legendre coefficients 8C(n). Solid line: Chebyshev coefficients 8C(n). 



If e is chosen as the backward error for a computed eigenvalue A, then the 
contour line of level e encloses all the complex numbers with the same backward 
error e for the pair (A, B). The larger the enclosed area, the worse-conditioned 
the eigenvalue. The diameter of the enclosed area gives an idea of the largest 
possible relative error on A. For a semi-simple eigenvalue, it is always possible to 
bound the error on A by the product of the condition number and the backward 
error. This is not possible for multiple defective eigenvalues and the spectral 
portrait is a useful alternative. 

On the example studied here, the computation is backward stable : we can 
then look at the contour line of level machine precision that is 10~ 16 (only 
"— log 10 e" appears on the figures). We see that this level curve encloses a large 
region of the spectrum, which tends to indicate a significant spectral instability 
in the matrix pairs under study. 

We display in figure [S] the spectral portrait for our eigenvalue problem using 
a resolution of L = 70 and N — 40 which corresponds to matrices of order 3150. 
We superpose the eigenvalues obtained using the QZ algorithm (black points) 
and the isolines of spectral portrait. For values of the spectral portrait larger 
than approximately 16 (lower part of the figure) the computed eigenvalues are 
completely undetermined in double precision. This corresponds to damping rates 
larger than 0.25 approximately. 

However, computation of pseudospectra is an expensive task and it is there- 
fore not feasible on production runs. There have been recent developments in 
the algorithms whereby one can obtain an approximation to the pseudospectra 
in a region near the interesting eigenvalues at reasonable cost J17I18J . However, 
those techniques must be used with special care as they are not totally reliable 
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Fig. 5. Spectral portrait. E = 1(T 4 , L = 70, N = 40. 
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Fig. 6. Plot of several eigenvalues ob- 
tained by perturbing randomly the two 
matrices A and B of eq. Q). The 
magnitude of the perturbation is the 
machine precision 2.22 x 10~ 16 . The 
shift is a fixed value near the the 
exact eigenvalue. The Ekman num- 
ber is E = 10 -4 and the resolution 
L = 94, N = 50. Each black dot 
in the plot is the difference between 
the computed eigenvalue and the ex- 
act one A = -0.3852109005 33277 + 
i0.65359 27894 40845. 
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Fig. 7. Plot of several eigenvalues ob- 
tained by making different calculations 
where the only change is the shift pa- 
rameter a of eq. ljl)l. a is changed by 
a random perturbation of magnitude 
10~ 5 near the exact eigenvalue. The 
parameters are the same as in figure 
El We remark that the two plots are 
hardly distinguishable. 
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Fig. 8. Same as figure^! only the calculations are done using extended precision. 



in the case where the matrix is highly nonnormal. A cheap technique that is 
used routinely to evaluate the sensitivity of eigenvalues to round-off error is to 
compute the eigenvalues of randomly perturbed matrices. This technique can 
be used without further coding on any eigenvalue solver, yet one must code the 
perturbations to the matrix elements. In the following we explore the impact of 
round-off errors by means of matrix perturbation, and present a new technique 
which gives the same results but does not require any coding at all. 

In figure [5] we plot the eigenvalues obtained by making several calculations 
on perturbed matrices. Each point in the figure is the eigenvalue obtained by 
perturbing the two matrices A and B of eq. J3J by random values uniformly 
distributed in the interval (— e m ,e m ), where e m — 2.22 x 10 -16 is the machine 
precision. The eigenvalues form a cloud of points concentrated in the neighbour- 
hood of the exact eigenvalue A = -0.38521 09005 33277 -M0.65359 27894 40845 
(this eigenvalue has been obtained using extended precision; it is the "exact" 
eigenvalue of the truncated problem © and not the one of the differential prob- 
lem (J2J). We did a statistical analysis on a large number of eigenvalues (50000). 
The results are summarized in tableland figure Both the real and imaginary 
parts follow quite well a Gaussian law: this can be seen in figureElwhere the prob- 
ability density functions of computed eigenvalues are plotted together with the 
gaussian curve that fits at best the data. More quantitatively, the skewness and 
kurtosis (table QJ are very close to those of the normal distribution (resp. and 
3). The order of magnitude of the round-off error is given by the standard devia- 
tion of the data which is a T ~ 7.69 x 10~ 6 for the real part and a^ ~ 6.83 x 10 -6 
for the imaginary part. Note that the covariance a TU , — Xa =1 (tj— r)(w.i — uj) (like 
the correlation coefficient p rul = (Xt " ) is small but non-zero, which means that 
the error distributions for the real and imaginary parts are slightly correlated. 
We remark that the standard deviations a T and a^ have similar values, i.e. the 
round-off error on r is of the same order of magnitude as that on to, even though 

\t\ « UJ. 

The standard deviations o r and a^ turn out to be essentially independent of 
the number of Chebyshev polynomials and spherical harmonics, provided that 
both spectra are well resolved. The values increase when the damping rate of the 



mode is increased, in perfect accordance with the plateaux observed in figures|21 
and|3J 





Matrix perturbation 


Shift perturbation 


Matrix perturbation in quad. prec. 


T 


-0.38520 966 


-0.38521089 


-0.3852109005 23 


To 


0.65359 249 


0.65359 283 


0.65359 2789414 


T— TQP 


1.24 x 10~ b 


8.67 x 10" a 


1.04 x 10" 11 


LU — LOQp 


-2.95 x 10 - ' 


3.67 x 10 _b 


-2.65 x 10~ XI 


Or 


7.69 x 10~ b 


7.80 x 10~ b 


4.99 x 10~ y 


Ou 


6.83 x 10~ b 


6.91 x 10~ b 


4.68 x 10~ y 


pTU> 


0.173 


0.142 


-0.476 


skewness(r) 


-0.016 


-0.004 


0.004 


skewness(tj) 


-0.003 


-0.011 


0.0004 


kurtosis(r) 


2.96 


2.98 


2.93 


kurtosis(tj) 


2.98 


2.98 


2.95 



Table 1. Statistics for the computed eigenvalues of figures and [7| We give 
the values of the averages r and to, standard deviations a T and <r w , cross corre- 
lation Ptui skewness and kurtosis for the perturbed matrix case (column 2), the 
perturbed shift case (column 3) and the perturbed matrix case using extended 
precision. 

(t qp = -0.38521 09005 33277, u QP = 0.65359 27894 40845) stands for the "exact" 
eigenvalue computed with quadruple precision 



In a second series of 50000 computations we did not perturb the matrices 
A and B but instead we perturbed the value of the Arnoldi-Chebyshev shift 
a by a small random quantity near the exact eigenvalue. With this method 
there is no need to modify the source code for the eigenvalue solver and/or for 
the construction of the matrices: we only need to change the shift parameter 
a on input to the eigenvalue solver. We obtain a cloud of eigenvalues (figure 
UJ| which looks almost identical to that obtained in figure Each point in the 
figure is the eigenvalue obtained by changing the real and imaginary part of 
the shift around the exact eigenvalue by random values uniformly distributed in 
the interval (— 10~ 5 , 10~ 5 ). Actually we have verified that the statistics does not 
depend on the amplitude of the shift perturbation. So there is no need to know 
a priori the exact value of the eigenvalue: any value of the shift which delivers 
the wanted eigenmode is good. The statistical values in tablenconfirm that the 
statistics of the eigenvalues obtained in the two approaches are almost the same. 

In a third series of 50000 computations the matrices A and B are perturbed 
as in the first series by random values uniformly distributed in the interval 
(— e m ,e m ), where e m = 2.22 x 10~ 16 . However in this series the computation is 
performed using extended precision. Thus we measure directly the sensitivity of 
the eigenvalues to perturbation of the matrices; in other words we compute the 
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Fig. 9. Probability density functions for the computed eigenvalues of figures HJ1 
and In abscissas are the differences between the real part of the eigenvalue 
and the average value r for the preturbed shift case (squares) and the perturbed 
matrix case (circles), and the differences between the imaginary part of the 
eigenvalue and the average value w for the preturbed shift case (stars) and the 
perturbed matrix case (plus). The continuous and broken lines corresponds to 
the gaussian curves which fita at best the data. We see that gaussian fit is almost 
perfect. 



e m -pseudoeigenvalue. From table Q we see that there are about three digits of 
difference between the standard deviations of the first two series and the present 
one: this means that the Arnoldi-Chebyshev algorithm adds an extra factor of 
order of magnitude 10 3 to the round-off error. 

6 Conclusions 

We have analyzed in this paper the errors that arise from the discretization and 
numerical computation of partial differential eigenvalue problems yielding large 
matrices. We have chosen as a model problem the two-dimensional eigenvalue 
problem yielded by the computation of inertial modes in a spherical shell. 

We have solved this problem using spectral methods for discretization and 
the incomplete Arnoldi-Chebyshev algorithm for solving the eigenvalue problem. 
The combination of these methods provides an efficient solver for these large 
(two-dimensional) eigenvalue problems. 

We have shown that the convergence of the eigenvalue and the eigenvector, 
with respect to spatial truncation, are tightly related: the absolute error of the 
eigenvalue decreases linearly with the Chebyshev truncation error, and quadrati- 
cally with the spherical harmonics truncation error, until round-off error becomes 
dominant. 

We found that most-damped modes are the most ill-conditioned and are 
therefore more sensitive to round-off error. This is made clear by the spectral 
portrait of the linear operator. Its computation is however very expensive and can 
be done only on small test problems. We have shown that a good estimation of 
the round-off error can be done in practice by performing different computations 
changing only the value of the Arnoldi-Chebyshev shift parameter on input; there 
is no need to do extra coding and/or to use external tools. It turns out that the 
round-off error on eigenvalues has an almost normal distribution, a result which 
can be used to reduce this kind of error. If the computation of a single eigenmode 
is not too expensive one could take advantage of this distribution of errors and 
perform N computations with random shifts; one can thus reduce the round-off 
error of the estimated eigenvalue by a factor yN. 
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